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Abstract. 

We present a joint implementation of dynamical-mcan-field theory (DMFT) with 
the pseudopotential plane- wave approach, via Wannier functions, for the determination 
of the electronic properties of strongly correlated materials. The scheme uses, as input 
for the DMFT calculations, a tight-binding Hamiltonian obtained from the plane- 
wave calculations by projecting onto atomic-centered symmetry-constrained Wannier 
functions for the correlated orbitals. We apply this scheme to two prototype systems: a 
paramagnetic correlated metal, SrV03, and a paramagnetic correlated system, V2O3, 
which exhibits a metal-insulator transition. Comparison with available Lincar-Muffin- 
Tin-Orbital (LMTO) plus DMFT calculations demonstrate the suitability of the joint 
DMFT pscudopotcntial-plane-wave approach to describe the electronic properties of 
strongly correlated materials. This opens the way to future developments using the 
pseudopotcntial-plane-wave DMFT approach to address also total-energy properties, 
such as structural properties. 
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1. Introduction 

Dynamical electron correlations play an important role in the physics of strongly 
correlated electronic materials, in particular in determining the properties of their 
paramagnetic and ferromagnetic phases, as well as their metal-insulator transitions. The 
latter transitions are related to some of the most dramatic effects observed in transition- 
metal oxides, such as the colossal magnetoresistance effect in doped manganites [H12]- 
The properties of transition-metal oxides are also known to be controlled by a strong 
and complex interplay between electronic, magnetic, and structural degrees of freedom 
[2]. This interplay leads to giant responses to small changes in external parameters 
such as temperature, pressure, magnetic field or doping, which makes such materials 
attractive for technological applications. 

A new theoretical framework has made possible in recent years the incorporation 
of dynamical correlations in electronic structure calculations of strongly correlated 
materials [3j EJ [5] . It combines the dynamical- mean- field theory (DMFT) of many-body 
physics with density-functional electronic structure calculations in the local-density 
approximation (LDA), or in the generalized gradient approximation (GGA), and is 
commonly referred to as LDA+DMFT. So far LDA+DMFT computations have been 
generally implemented with the Linear Muffin-Tin Orbital (LMTO) method, and in 
some cases also with mixed-basis methods [6], which incorporate also atomic-like orbitals 
in the basis sets. In view of the importance, however, of the interplay between structural 
and electronic properties in correlated oxides, it would be highly desirable to have a 
joint implementation of DMFT with the electronic-structure pseudopotential plane- 
wave method. The latter method is indeed known to be well suited to address also 
total-energy properties, such as structural properties. This is related to the plane-wave 
basis, in particular, which does not depend on the atomic positions in the unit cell and 
whose completeness can be controlled by a single parameter, the plane-wave kinetic 
energy cutoff. 

Recently, approaches based on Wannier functions have been proposed to carry 
out LDA+DMFT computations, using e.g., atomic-centered Wannier functions [7], or 
Maximally localized Wannier functions [6], or Wannier functions constructed with the 
N-th order muffin-tin-orbital method [Sj. For strongly correlated materials, Wannier 
functions represent a convenient, physically sound set of localized orbitals for the 
correlated electrons that can be used to construct an interface between DMFT and 
LDA/GGA calculations. This is especially true for LDA/GGA methods which do not 
employ atomic-like basis functions, as is the case of the pseudopotential plane-wave 
method. 

Here we present a joint implementation of DMFT with the pseudopotential plane- 
wave approach, and demonstrate its suitability to determine the electronic properties 
of correlated oxides. This is a prerequisite to address, in the future, also total- 
energy and structural properties. We use as input for the DMFT calculations a 
tight-binding Hamiltonian, B. TB (k), constructed from the pseudopotential plane- wave 
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calculations using atomic-centered symmetry-constrained Wannier functions for the 
correlated orbitals. We apply this approach to two different test cases: a correlated 
paramagnetic metal, SrVC>3, with a simple cubic perovskite structure, and a correlated 
paramagnetic system, V2O3, which has a more complex trigonal crystal structure, and 
exhibits a Mott-Hubbard transition. To our knowledge, this is the first successful 
implementation of LDA+DMFT using the pseudopotential plane-wave approach. 

2. Method 

2.1. LDA+DMFT calculation scheme 

The standard (although simplified [9]) scheme generally used to carry out LDA+DMFT 
calculations includes the following three steps. 

First, an LDA/GGA self-consistent electronic-structure calculation is performed 
for the crystal and phase of interest. The calculated bands associated with the relevant 
correlated orbitals (e.g., the transition- metal 3d-t2 9 orbitals in our examples) are then 
mapped onto a k-dependent tight-binding Hamiltonian, Hl^ A (k), where k is a vector 
of the crystal Brillouin zone (BZ). In practice, in the case of LDA/GGA calculations 
performed using, as a basis set, atomic-like orbitals, this Hamiltonian is readily obtained 
as: 

HlSaO*) — E ^mm' ' ) (k)cj cmtT Ckm'cr (1) 

m,m',<T 

where H^jf^ (k) are the matrix elements, in a given k and a spin-polarization subspace, 
of the LDA/GGA Hamiltonian between the relevant correlated atomic-like orbitals, 
indexed by m, m', and cl^ mo . (ck mcr ) is the creation (annihilation) operator for an electron 
in orbital m with spin o and wavevector k. As we are interested here in the case of 
paramagnetic phases, we will omit in the following the spin label a in H^f\ We 
note that, at the LDA/GGA level, the paramagnetic phase is simply modeled by a 
non- magnetic state. We are including in the sum, in Eq. (JTJ), a minimal basis set of 
m, m' orbitals corresponding to the relevant onsite correlated orbitals near the Fermi 
energy. More generally, however, additional non-correlated orbitals may be included, 
when needed, to improve the description of the spectrum, e.g., further away from the 
Fermi energy [10] . For simplicity, we will focus in the following on the case of periodic 
systems with a single type of correlated atomic sites. 

Next, the LDA tight-binding Hamiltonian is supplemented with onsite Coulomb 
interactions for the correlated orbitals, in a many-body Hamiltonian of the form: 

k,m,m' ,cr I,m,cr 

+ ^Un Im1 n Iml + E (V-L'J)n; mff n/mV- (2) 

I,m I ,m^m' ,a,a' 

The index I labels the lattice sites, m, m! label the relevant correlated atomic orbitals, 
and cr, a' the spin states (f and |), and n/ mCT is the operator for the occupation of the 
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correlated orbitals. The first term, on the right-hand side of Eq. (J2J), is the LDA part of 
the Hamiltonian. The second term is a double- counting correction, formally introduced 
to remove the onsite Coulomb interactions already present, in an average way, in the 
LDA Hamiltonian. In principle, the double-counting potential Ae dc may be taken as 
[I]: Ae dc = U{n - §), where n is the average occupation per correlated orbital. In 
practice, however, if one takes into account only d orbitals (or only / orbitals), the 
double-counting correction acts on the whole d (/) band and shifts it by Ae dc . This 
potential amounts thus to a rigid energy shift of the quasiparticle spectrum in the DMFT 
calculations. It can therefore be absorbed into the chemical potential, which in turn is 
determined by the number of electrons. Hence the actual value of Ae dc has no influence 
on the spectral properties. The third and forth terms, on the right-hand side of Eq.([2]), 
are the interaction terms, where U is the onsite Coulomb repulsion parameter, J the 
Hund's rule exchange parameter, and V = U — 2 J [TT] . 

Finally, the model Hamiltonian, in Eq.([2]), is solved by means of DMFT. DMFT 
maps the lattice model onto an effective single-impurity problem subject to a self- 
consistent condition on the impurity self-energy, £(e), or, equivalently, on the local 
Green's function, G(e) [12] • This mapping represents an exact solution in the limit 
of infinite dimension of the lattice problem. The local impurity Green's function and 
self-energy matrices are related through: 

<w(e) = 7^—1 dk([(e - - H TB (k) - S(e)]~ 1 )w, (3) 

11 BZ J 

where m,m' label the correlated orbitals, \i is the chemical potential, 1, S(e), and 
H TjB (k) are the unitary, self-energy, and LDA tight-binding-Hamiltonian matrices, 
respectively, and the integration extends over the BZ with volume Qbz- Several different 
approaches may be used to solve the effective impurity problem, including numerical 
renormalization group, exact diagonalization, noncrossing approximation, and Quantum 
Monte Carlo (QMC). In the present work, we employ, as impurity solver, the auxiliary 
fields QMC method by Hirsch and Fye [13] . For a given input impurity self-energy and 
Green's function matrices, the impurity solver yields a new set of Green's function and 
self-energy matrices [12] . The DMFT equations are solved in an iterative self-consistent 
cycle, until self-consistency is reached. 

2.2. Tight-binding Hamiltonian from pseudopotential plane-wave calculations 

Wannier functions provide a physically sound basis set to construct model Hamiltonians 
for correlated electrons. In the case of LDA/GGA plane- wave calculations, they also 
represent a practical route to build Hl^ A (k), given the delocalized nature of the basis 
functions. 

The Wannier functions are obtained here using the approach by Ku et al. [T4"]- 
inspired from the work by Marzari and Vanderbilt [T5] . to obtain atomic-centered 
Wannier functions with a given symmetry. First, a set of M trial functions are 
generated from the (pseudo) atomic wave-functions of the correlated orbitals: |y3 m k) = 
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X)r e* kR |<^), where R are the lattice translation vectors and |(/? R ) are the atomic 
(pseudo) wave-functions, in the unit cell R, with a given symmetry. A set of M non- 
orthogonal Wannier functions (WF) in k space (or tight-binding Wannier functions) is 
then obtained by projecting these trial functions onto a set of Bloch functions, l^jk), 
belonging to a chosen (correlated bands) subspace: 

IWmk) = (^ikl^mOl^k)- ( 4 ) 

j=Ni 

The sum is over a Bloch subspace defined by imposing either some fixed band 
numbers, N± < j < N2, or an energy window, E\ < ej(k) < E2, for the electronic 
bands e 3 -(k). These Wannier functions are then orthogonalized, by diagonalizing 
their overlap matrix O mm >, yielding a set of M orthogonalized Wannier functions: 
|H / mk ) = E m -(0" 1/2 ) mm # m ' k ). 

Using this basis set of Wannier functions, the tight-binding Hamiltonian matrix is 
given by: 

Hll> WF W = E (W mk |^ k )fe|M/ m , k )6 j (k), (5) 

3=Ni 

and can be used, in place of i?^,(k), in the many-body Hamiltonian, in Eq. (J2J). 
It should be noted that the conventional, real-space (lattice-site related) Wannier 
functions, (H 7 ^), are simply the Fourier transforms of the tight-binding Wannier 
functions: \W^) = j^- J2k e~ ikR |W / mk)- Only the tight-binding Wannier functions, 
however, are explicitly needed in the implementation. 

We implemented this scheme to construct H TjB (k), in k space, in the framework 
of the pseudopotential plane-wave method. This was done as an interface between 
the Quantum-Espresso pwscf package [16] and the QMC-DMFT code. The matrix 
elements (^kl^mk), in Eq. (BJ, can be conveniently evaluated in reciprocal space, in the 
pseudopotential plane- wave scheme. For the trial wave-functions |y mk ), we generated 
wave-functions belonging to the point-group representations of the correlated atomic 
site. This was done by diagonalizing the occupation matrix [17] — calculated in an 
initial (arbitrary) orthogonal basis set in the I angular-momentum subspace of the 
correlated atomic orbitals, and using the corresponding eigenstates, which belong each 
to a specific representation. The Hamiltonian H TjB (k) was evaluated on a k-point grid 
in the irreducible part of the BZ. The integration, in Eq. (jHJ), was performed using the 
analytical tetrahedron method [IE], and restricted to the irreducible part of the BZ by 
symmetrizing the Green's function matrix. 

2.3. Technical details: QMC-DMFT computations 

Computationally, the most involving part of the calculations is the evaluation of the 
path integral in the auxiliary-field QMC method, to solve for the local impurity Green's 
function [12]. The QMC method maps the interacting electron problem onto a sum 
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of non-interacting problems by means of the Trotter discretization and Hubbard- 
Stratonovich transformation [T21 EES]. The imaginary time integrals are represented, 
in the Trotter discretization, by L imaginary-time slices of size A = (3/L, with f3 = 
1/KgT. For an M-orbital impurity problem, the Hubbard- Stratonovich transformation 
introduces M(2M — 1) auxiliary Ising fields for each time slice. In addition to the 
tolerance parameter of the DMFT self-consistency, the number of time slices L and the 
number of Monte Carlo sweeps Nmc f° r the stochastic integration of the path integral 
are the sole convergence parameters of the QMC-DMFT calculations. It should be noted 
that the computational cost of the QMC algorithm scales, to the leading order in L, as 
~ M(2M - 1)N MC L 3 p]. 

The QMC method has the advantage of having formally no approximation. Very 
low temperatures, however, are not accessible, as the numerical effort scales as 1/T 3 . 
From the imaginary-time self-consistent Green's function obtained from the QMC- 
DMFT computations, the real-frequency single particle spectral functions are computed 
using the maximum entropy method [ZD] . 

3. Applications 

The joint DMFT pseudopotential plane-wave scheme described in the previous sections, 
was applied to two test cases, SrV0 3 and V 2 3 . For the density functional calculations, 
we used the Perdew-Burke-Ernzerhof exchange-correlation functional [21] together with 
Vanderbilt utrasoft pseudopotentials [22] . We used a kinetic energy cutoff of 35 Ry (350 
Ry) for the plane- wave expansion of the electronic states (core-augmentation charge). 
The self-consistent calculations were performed with a (4,4,4) Monkhorst-Pack k-point 
grid [23]. For the computation of H TB (k) and of the Green's functions matrix, in Eq. 
([3]), we used a (10,10,10) k-point grid centered at T. The experimental values of the 
lattice parameters of SrV0 3 (a = 3.84 A) [23] and V 2 3 (a = 4.95 A, c = 14.00 A) [25] 
have been used in our calculations. 

For SrV0 3 , we set the onsite Coulomb interaction U = 5.55 eV and Hund's rule 
parameter J = 1 eV [26]. For V20 3 , we used J = 0.93 eV [25] and several different values 
of U 0. The QMC-DMFT calculations were performed at T = 580 K (j3 = 20 eV" 1 ), 
using 80 imaginary-time slices. In the case of V20 3 , we also performed calculations 
at T = 1160 K (P = 10 eV" 1 ), using 40 imaginary-time slices. In all QMC-DMFT 
calculations we used ~ 10 6 Monte Carlo sweeps. 

3.1. SrV0 3 

SrV0 3 is a prototype d 1 correlated paramagnetic metal. It has a simple cubic perovskite 
structure, and remain paramagnetic down to low temperatures. It is an ideal test case 
for first-principles many-body calculations. In Fig. [IJ we show the density of states 
(DOS) obtained for SrV0 3 from the GGA pseudopotential calculations. The spectrum 
is in agreement with previous LDA calculations [HI EE]. The valence states of SrV0 3 
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consist of completely occupied oxygen 2p states, located in the energy range -7 eV to 
-2 eV below the Fermi energy, and partially occupied V-3d t<i g states, near the Fermi 
energy. 

From the pseudopotential plane-wave calculations, we generated the V-3d- t 2g and 
e g Wannier functions from the Bloch functions corresponding to the 5 lowest-energy (3d) 
bands, in the energy window -1 to 5.5 eV. The corresponding projected Wannier t 2g and 
e g DOS's are also displayed in Fig. [TJ Because of the ideal octahedral symmetry of the 
V sites, hybridization is forbidden between the t2 g and e g states. H TjB (k) is hence block 
diagonal with respect to these two subspaces. We have used, as correlated subspace, 
the t2 g Wannier-functions subspace, and the corresponding t2g-Hamiltonian block for 
the DMFT computations. 



2p 



V3d 



10 



8 



-- V3dt 
--V3de, 



2g 




Figure 1. (Color online) Density of states of SrVC>3 obtained from the GGA 
pseudopotential plane-wave calculations. The short- (long-) dashed line shows the V- 
3d- tig (e g ) Wannier-projected density of states. The generic characters of the bands 
are also indicated at the top of the figure. The zero of energy corresponds to the Fermi 
level. 



It should be noted that, in the special case of cubic symmetry, the Green's function 



matrix of the t 2g states, in Eq. ([3]), may be expressed as: G mm '(e) = J 



deD{e') 



<5n 
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This is valid, however, only when the local matrices (self-energy and Green's function 
matrices) are proportional to the unitary matrix, i.e., in the case of cubic and 
higher symmetry. In the present work we always used the more general Hamiltonian 
formulation with the k space integration. 

In Fig. El we display the corresponding single-particle spectral function obtained 
from the QMC-DMFT calculations at T = 580 K. Taking into account the correlations 
effects within the t2 9 manifold leads to substantial modifications in the single-particle 
spectrum relative to the GGA result. Correlations effects are responsible for a 
lower Hubbard band around -2 eV, an upper Hubbard around 2.5 eV, and a well- 
pronounced quasiparticle peak at the Fermi energy. This is in general agreement with 
the photoemission and inverse photoemission experiments on SrV03 [27] . The results in 
Fig. El compare well with previous LMTO-based LDA+DMFT computations performed 
with the same value of U [26], and are also consistent with the mixed-basis LDA+DMFT 
calculations using somewhat smaller values of U [6]. 

3.2. V 2 3 

V2O3 is a Vanadium d 2 system. The high-temperature paramagnetic phase of V2O3 has 
a corundum trigonal crystal structure, with 4 equivalent V sites in the unit cell. Within 
the corundum structure, each V ion is surrounded by a distorted oxygen octahedron 
[28| [29] , The V ions are arranged in pairs along the c axis, with a stacking that can 
be obtained, starting from an ideal chain of V ions along c, by introducing vacancies 
at every third site [25]. While the V-V pairs along the c axis are surrounded by face- 
sharing oxygen octahedra, in the a-b plane, each V ion has three nearest neighbors with 
edge-sharing oxygen octahedra [28] . 

Assuming ideal V octahedral sites in V2O3, the V-3d atomic states are split into t2 g 
and e g orbitals, where the two degenerate e g orbitals are empty and the three degenerate 
t2 g orbitals are filled with 2 electrons. When the symmetry is further reduced by a small 
trigonal distortion of the octahedra in the corundum structure, the t2 9 orbitals further 
splits into a non-degenerate a± g orbital oriented along the c axis, and two degenerate 
orbitals oriented predominantly in the a-b plane. 

In Fig. El we show the DOS of V2O3 obtained from the GGA pseudopotential 
calculations. The oxygen 2p states of V2O3 are located roughly between -8 and -4 
eV below the Fermi energy. The V-3d ^Tike states are located in an energy window 
between -1.5 eV and 1.5 eV, and the e a g states are approximatively between ~ 2 eV and 
4 eV. This is consistent with previous LDA calculations [25j EH ED]. The GGA (and 
also LDA) calculations yield a metallic phase for V2O3, with a high DOS at the Fermi 
energy. Experimentally, instead, the paramagnetic corundum phase is found to be an 
insulator at low pressure. 

We generated Wannier functions with a\ g and symmetry from the Bloch states 
enclosed in the energy window -1.5 to 1.5 eV (see Fig. [3]). The corresponding Wannier 
projected DOS's are displayed in Fig. [3j We used these Wannier functions to construct 
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E(eV) 

Figure 2. (Color online) Spectral function for SrVC>3 at T = 580 K obtained 
from the GGA+DMFT computations using the ti g projected Hamiltonian from the 
pscudopotcntial plane- waves calculations. The zero of energy corresponds to the Fermi 
level. 

H TjB (k) and carry out the GGA+DMFT computations. The results are shown in Figs. 
H] and [5] for two different temperatures: T = 1 160 K and 580 K, respectively. The 
separate a\ g and contributions to the spectral function are displayed in the upper 
and lower panels, respectively. At T = 1160 K, calculations were performed for U = 5.6 
and 6 eV. An insulating phase is obtained for U = 6 eV, while the system is still 
metallic at U = 5.6 eV. The results in Fig. H] agree well with the available LMTO- 
based LDA+DMFT calculations performed at the same temperature [25] [30] , and in 
particular with the results obtained using a Wannier-projected Hamiltonian [7]. The 
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2p V3dt 2 V3de 




E(eV) 

Figure 3. (Color online) Density of states of V2O3 obtained from the GGA 
pseudopocntial plane-waves calculations. The short- (long-) dashed line shows the 
V-3d- (dig) Wannier-projected density of states. The generic characters of the 
bands are also indicated at the top of the figure. The zero of energy corresponds to 
the Fermi level. 

pseudopotential implementation is found to yield a slightly larger critical U for the 
insulating phase (U ~ 6 eV) compared to the LMTO implementation (U ~ 5.5 eV) 
[7j. This is related to the band width of the LDA/GGA i2gTike states, which is slightly 
larger in the pseudopotential case compared to the LMTO case. 

At T = 580 K, we performed calculations for U = 5 and 5.6 eV. An insulating phase 
is found at U ~ 5.6 eV, whereas at U = 5 eV the system is metallic, with a large DOS 
at the Fermi energy. In the latter case, one observes a quasi-particle peak at the Fermi 
energy in the a\ g orbital-resolved spectral function. The results at T = 580 K, in Fig. 
are in good qualitative agreement with the results of LDA+DMFT calculations using 
a Wannier Hamiltonian constructed with the N-th order muffin-tin-orbital method [31 J . 
In the latter study, a critical U of 4.2 eV was found for the insulating phase. The larger 
critical U obtained here (5.6 eV) is attributed mainly to the difference in the crystal 
structure; in Ref. [31] the (Vo.962Cr ao38 ) 2 03 crystal structure was considered — which 
corresponds to the experimental insulating structure. The increased U is also due in 
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part to the broader GGA band found in the pseudopotential calculations. One observes 
also some small differences in the peak structure, between the present results in Fig. 
[5] and the spectral function in Ref . [31] . These are attributed mainly to differences in 
the DMFT calculational details, such as the use of two different interpolation schemes 
(Ulmke-Janis-Vollhardt scheme [32] versus cubic splines) for the Fourier transformation 
of the local Green's function, in the two studies. 
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Figure 4. (Color online) V-3<i ai s and contributions to the spectral function of 
V 2 3 at T = 1160 K, as obtained from the GGA+DMFT computations using the t 2g - 
like Hamiltonian constructed from the pseudopotential plane-waves calculations. The 
results are shown for two different values of the onsitc Coulomb interaction, U = 5.6 
and 6 eV. The vertical dotted line indicates the Fermi energy. 



3.3. Discussion and outlook 

The results we have presented here show that the DMFT plus pseudopotential plane- 
wave scheme is both a practical and suitable approach for the determination of the 
electronic properties of strongly correlated oxides. A promising extension of this 
approach concerns the determination of total-energy properties, and in particular of 
structural properties of correlated systems. The total energy, within the LDA+DMFT, 
can be expressed as [33]: E = E LDA - E DC - E m ,k e m^k A + ( h lda) + ( H u), where E LDA 
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is the LDA total energy, E DC is the double-counting energy corresponding to the second 
term on the right-hand side of Eq. (j2J), Em.kfi' 1 i s the sum of the LDA valence-state 
eigenvalues, (Hl^ A ) = tr[H TB G], and (Hjj) is the interaction energy, corresponding 
to the third and fourth term on the right-hand side of Eq. (j2J), computed from the 
double occupancy matrix [34]. The application of the DMFT pseudopotential plane- 
wave approach to the determination of structural relaxations is a line of development 
we are currently pursuing. 



1 1 I 1 I 1 I r 




E (eV) E (eV) 

Figure 5. (Color online) Same data as in Fig. 4, but for T = 580 K and U = 5 and 
5.6 eV. 

4. Summary and conclusions 

We have presented an implementation of the LDA+DMFT approach within the 
pseudopotential plane-wave framework. This scheme was applied to two different 
test cases, SrVC>3 and V2O3. Comparison with available LMTO-based LDA+DMFT 
calculations demonstrated the suitability of the joint DMFT pseudopotential-plane-wave 
scheme to describe the electronic properties of strongly correlated materials. This opens 
the way to future developments using the DMFT pseudopotential-plane-wave approach 
to address also total-energy and hence structural properties of correlated systems. 
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